




library(mvtnorm)
library(arm)
library(BRugs)
library(R2OpenBUGS)
library(coda)
library(car)
library(foreign)
library(vioplot)



setwd("~/../Dropbox/iREDS/Data/Analysis/scales/")



# Read in OpenBUGS results and create figure

setwd("~/../Dropbox/iREDS/Data/Analysis/scales/S1/")
coda.data<-read.openbugs(stem="")
results.S1<-as.data.frame(rbind(coda.data[[1]], coda.data[[2]], coda.data[[3]]))

setwd("~/../Dropbox/iREDS/Data/Analysis/scales/S2/")
coda.data<-read.openbugs(stem="")
results.S2<-as.data.frame(rbind(coda.data[[1]], coda.data[[2]], coda.data[[3]]))

setwd("~/../Dropbox/iREDS/Data/Analysis/scales/S3/")
coda.data<-read.openbugs(stem="")
results.S3<-as.data.frame(rbind(coda.data[[1]], coda.data[[2]], coda.data[[3]]))

setwd("~/../Dropbox/iREDS/Data/Analysis/scales/S5/")
coda.data<-read.openbugs(stem="")
results.S5<-as.data.frame(rbind(coda.data[[1]], coda.data[[2]], coda.data[[3]]))

setwd("~/../Dropbox/iREDS/Data/Analysis/scales/S6/")
coda.data<-read.openbugs(stem="")
results.S6<-as.data.frame(rbind(coda.data[[1]], coda.data[[2]], coda.data[[3]]))

setwd("~/../Dropbox/iREDS/Data/Analysis/scales/S8/")
coda.data<-read.openbugs(stem="")
results.S8<-as.data.frame(rbind(coda.data[[1]], coda.data[[2]], coda.data[[3]]))

coda.data.merge<-as.mcmc(rbind(coda.data[[1]], coda.data[[2]], coda.data[[3]]))
sink("HPDintervals.txt")
HPDinterval(coda.data.merge)
sink()





setwd("~/../Dropbox/iREDS/Data/Analysis/scales/")


# extract treatment effect estimates and standardize for Cohen's D
std.deltatheta<-c(0.86, 1.17, 1.33, 1.06, 1.16, 1.17)
alpha<- as.data.frame(cbind(results.S1$S1.alpha1/std.deltatheta[1], results.S2$S2.alpha1/std.deltatheta[2], results.S3$S3.alpha1/std.deltatheta[3], 
              results.S5$S5.alpha1/std.deltatheta[4], results.S6$S6.alpha1/std.deltatheta[5], results.S8$S8.alpha1/std.deltatheta[6]))
names(alpha)<-c("S1.alpha1", "S2.alpha1", "S3.alpha1", "S5.alpha1", "S6.alpha1", "S8.alpha1")


## Main figure

nrows<-length(alpha[,1])
ncols<-6 # this is the number of rows in the figure
figdata<-as.data.frame(matrix(rep(NA, nrows*ncols*2), nrow = nrows*ncols, ncol=2))

for (j in 1:ncols) {
  for (i in 1:nrows) {
    if (j==1) {
      figdata[i+5*nrows,1] <- alpha$S8.alpha1[i]
      figdata[i+5*nrows,2] <- j}
    if (j==2) {
      figdata[i+4*nrows,1] <- alpha$S3.alpha1[i]
      figdata[i+4*nrows,2] <- j}
    if (j==3) {
      figdata[i+3*nrows,1] <- alpha$S2.alpha1[i]
      figdata[i+3*nrows,2] <- j}
    if (j==4) {
      figdata[i+2*nrows,1] <- alpha$S5.alpha1[i]
      figdata[i+2*nrows,2] <- j}
    if (j==5) {
      figdata[i+1*nrows,1] <- alpha$S6.alpha1[i]
      figdata[i+1*nrows,2] <- j}
    if (j==6) {
      figdata[i,1] <- alpha$S1.alpha1[i]
      figdata[i,2] <- j}
  }
}

figdata[,2] <- factor(
  figdata[,2],
  levels = c("1", "2", "3", "4", "5", "6"),
  labels = c("Preserve Replication Materials", "Reasons for Data Management Policy", "Reasons for Authorship Policy",
             "Lab Disagreement", "Respectful Discussion", "Relevance of Ethics Discourse")
)




library(ggplot2)
library(ggridges)
library(ggpubr)
library(viridis)

theme_set(theme_ridges())



## USE THIS ONE -- THIS IS THE VERSION IN THE PAPER!
ggplot(figdata, aes(x=figdata[,1], y=figdata[,2], fill=0.5 - abs(0.5-..ecdf..))) +
  stat_density_ridges(geom = "density_ridges_gradient", scale=0.9, calc_ecdf = TRUE) +
  scale_fill_viridis(name = "Prob(Tail)", direction = -1) +
  labs(
    title = 'Treatment Effect Estimates',
    subtitle = 'Cohen\'s D',
    x = "Posterior Distribution"
  ) +
  theme_ridges(font_size = 13, grid = TRUE) + theme(axis.title.y = element_blank())







# old violin plot -- we don't use this!!
op <- par(no.readonly = TRUE) # the whole list of settable par's.
par(mar=c(5,15,5,1), mfcol=c(1,1), xpd=FALSE, las=1)

# Draw the plot
with(alpha , vioplot( S8.alpha1, S3.alpha1, S2.alpha1, S5.alpha1, S6.alpha1, S1.alpha1,   col=topo.colors(6), 
                        range = 1.75, lwd=2, border = NA, lineCol = "gray",
                        names=c("Preserve Replication Materials", "Reasons for Data Management Policy", "Reasons for Authorship Policy", 
                                "Lab Disagreement", "Respectful Discussion", "Relevance of Ethics Discourse"), 
                        main=paste("Treatment Effect Estimates (Cohen's D)"), xlab=paste("Posterior Distribution"),
                        horizontal = TRUE))  
abline(v=0)
abline(h=3.5, lty=2)
text(x=-1.75, y=5.7, "Lab Communication")
text(x=-1.75, y=5.4, "Scales")
text(x=-1.75, y=2.7, "Training Content")
text(x=-1.75, y=2.4, "Scales")

par(op)





# rho plot

setwd("~/../Dropbox/iREDS/Data/Analysis/single_items")
coda.data<-read.openbugs(stem="")
results<-as.data.frame(rbind(coda.data[[1]], coda.data[[2]], coda.data[[3]]))
setwd("~/../Dropbox/iREDS/Data/Analysis/scales/")


rho<- as.data.frame(cbind(sample(results$S1.rho, length(results.S1$S1.rho)),
                          sample(results$S2.rho.1, length(results.S1$S1.rho)),
                          sample(results$S2.rho.2, length(results.S1$S1.rho)),
                          results.S1$S1.rho, results.S2$S2.rho, results.S3$S3.rho, 
                          results.S5$S5.rho, results.S6$S6.rho, results.S8$S8.rho))
names(rho)<-c("S.1.rho", "S.2.rho.1", "S.2.rho.2", "S1.rho", "S2.rho", "S3.rho", "S5.rho", "S6.rho", "S8.rho")


op <- par(no.readonly = TRUE) # the whole list of settable par's.
par(mar=c(5,15,5,1), mfcol=c(1,1), xpd=FALSE, las=1)

# Draw the plot
with(rho , vioplot( S8.rho, S3.rho, S2.rho, S5.rho, S6.rho, S1.rho, S.1.rho, S.2.rho.2, S.2.rho.1,
                    col=topo.colors(9), 
                      range = 1.75, lwd=2, border = NA, lineCol = "gray",
                      names=c("Preserve Replication Materials", "Reasons for Data Management Policy", "Reasons for Authorship Policy", 
                              "Lab Disagreement", "Respectful Discussion", "Relevance of Ethics Discourse",
                              "Discussion Changed Ethics Views", "Lab Has Data Management Plan", "Lab Has Authorship Policy"), 
                      main=paste("Posteriors for Rho"), xlab=paste("Posterior Distribution"),
                      horizontal = TRUE))  
abline(v=1)
abline(h=6.5, lty=2)
text(x=-1.75, y=5.7, "Lab Communication")
text(x=-1.75, y=5.4, "Scales")
text(x=-1.75, y=2.7, "Training Content")
text(x=-1.75, y=2.4, "Scales")

par(op)


